# Prelims #
rm(list = ls())

# Prelims #
library(dplyr)
library(sandwich)
library(MASS)

#load the data
d <- read.csv("pncc_dy_16042022.csv")
d <- d[d$Negotiation != 2, ] # remove implementation phase

# Figure 2: Rate of negotiation by type and region. ####
x <- d %>% dplyr::select(DyadId, Year, Negotiation, MedType, Region) %>% unique()

# bilateral talks
x$bi <- 0
x$bi[d$MedType == "Bilateral Talks"] <- 1

# un mediation
x$un <- 0
x$un[d$MedType == "UN Mediation"] <- 1

# non un
x$nonun <- ifelse(x$MedType == "Non-UN Mediation", 1, 0)

x <- x %>% group_by(Region) %>%  mutate(all = n()) %>% mutate(all_neg = sum(Negotiation)) %>% 
  mutate(bi_neg = sum(bi)) %>% mutate(med_un = sum(un)) %>% 
  mutate(med_nonun = sum(nonun)) %>%   ungroup()

x <- x %>% dplyr::select(Region, all, all_neg, bi_neg, med_un, med_nonun) %>% unique()
x$rat_allneg <- x$all_neg / x$all
x$rat_bineg  <- x$bi_neg / x$all
x$rat_medun  <- x$med_un / x$all
x$rat_mednonun  <- x$med_nonun / x$all

## order and labels
x <- x[order(x$rat_allneg, decreasing = T),]
x$Region[x$Region == 2] <- "Mid. East"
x$Region[x$Region == 1] <- "Europe"
x$Region[x$Region == 3] <- "Asia/Pacific"
x$Region[x$Region == 4] <- "Africa"
x$Region[x$Region == 5] <- "Americas"

a <- x %>% dplyr::select(rat_medun, rat_mednonun, rat_bineg,) %>% as.matrix()
a <- t(a)

### colors ####
percent <- 20
fark <- 0

# un #
col   <- "dodgerblue3"
col1  <- rgb(red = col2rgb(col)[1] - fark, 
             green = col2rgb(col)[2] - fark/2,
             blue = col2rgb(col)[3], max = 255,
             alpha = (100-percent)*255/100)

# non un #
col   <- "salmon"
col2  <- rgb(red = col2rgb(col)[1] - fark, 
             green = col2rgb(col)[2] - fark/2,
             blue = col2rgb(col)[3], max = 255,
             alpha = (100-percent)*255/100)


# bilateral#
col   <- "grey87"
col3  <- rgb(red = col2rgb(col)[1] - fark, 
             green = col2rgb(col)[2] - fark/2,
             blue = col2rgb(col)[3], max = 255,
             alpha = (100-percent)*255/100)

dev.off()
png("regions.png", units="cm", 
    width = 12, height = 9, res = 600)
par(mfrow = c(1,1), mar = c(1.25,1.25,1.25,0), oma = c(0, 0, 0, 0), mgp = c(0,0.25,0), cex.axis = 0.75, tck = -0.01)
barplot(a, names.arg = x$Region, las = 1, ylim = c(0,0.5), col = c(col1, col2, col3), 
        main = "Rate of Negotiation", 
        cex.main = 0.75, cex.names = 0.6, cex.lab = 0.5)
legend("topright", legend=c("UN Mediation", "Non-UN Mediation", "Bilateral"), pch = c(15,15,15), col=c(col1, col2, col3), bty = "n",
       pt.cex = 2, cex = 0.75)
dev.off()

# Figure 3: Time trends in civil conflict negotiations.####
x <- d %>% dplyr::select(DyadId, Year, Negotiation, Mediation) %>% unique()

x <- x %>% group_by(Year) %>% mutate(all = n()) %>% mutate(neg_sum = sum(Negotiation)) %>% mutate(med_sum = sum(Mediation)) %>% ungroup()
x <- x %>% dplyr::select(Year, all, neg_sum, med_sum) %>% unique()
x <- x[order(x$Year),]

#colors
col1   <- "grey87"
colmed <- "salmon"
colneg <- "dodgerblue3"

##rate
x$rate <- x$neg_sum/ x$all
x$rate_neg <- (x$neg_sum - x$med_sum)/ x$all
x$rate_med <- x$med_sum/ x$all

length(unique(d$DyadId))
length(unique(d$GWNoA))

table(d$Exit)
sum(d$OtherDocument)

#### plot trends##
dev.off()
png("trends.png", units="cm", 
    width = 16, height = 9, res = 600)
par(mfrow = c(1,1), mar = c(1,1,0,0), oma = c(0, 0, 0, 0), mgp = c(1,0.2,0), cex.axis = 0.75) 
plot(x$Year, x$all, ylim = c(0, 80), t = "n", axes = F, xlab = "", ylab = "")
axis(2, at = seq(0, 80, 10), las = 1, tck = -0.01)
par(mgp = c(1,0.075,0), cex.axis = 0.75) 
axis(1, at = c(seq(1975, 2010, 5), 2013), tck = -0.01)
polygon(x = c(x$Year, rev(x$Year)), y = c(rep(0, nrow(x)), rev(x$all)), col = col1, border = NA)
polygon(x = c(x$Year, rev(x$Year)), y = c(rep(0, nrow(x)), rev(x$neg_sum)), col = colneg, border = NA)
polygon(x = c(x$Year, rev(x$Year)), y = c(rep(0, nrow(x)), rev(x$med_sum)), col = colmed, border = NA)
#lines(x$Year, x$all)
#lines(x$Year, x$med_sum)
#lines(x$Year, x$neg_sum)
legend("topleft", legend=c("All Dyads", "All Negotiation", "Only Mediated"), pch = c(15,15,15), col=c(col1, colneg, colmed), bty = "n",
       pt.cex = 2, xjust = 1, adj = c(0,0.5), cex = 0.8)
dev.off()

# Figure 6: Rebel strength and negotiation type. Estimates based on multinomial logit (Model 2). ####
## rebel strength
d <- read.csv("rebstr.csv")

b <- d[d$type == "negotiation",]
m <- d[d$type == "mediation",]

##plot
dev.off()
png("rebstr.png", units="cm", 
    width = 8, height = 10, res = 600)
par(mfrow = c(1,1), mar = c(1.6,2,1,1.6), oma = c(0, 0, 0, 0), mgp = c(2,0.25,0), cex.lab = 0.7, cex.main = 0.9)
plot(1:3, m$pe, pch = 15, ylim = c(0,0.4), axes = F, xlab = "", ylab = "", col = "salmon",
     main = "Rebel Strength")
lines(1:3, m$pe, col = "salmon")
segments(x0 = 1:3, x1 = 1:3, y0 = m$lb, y1 = m$ub, col = "salmon")
points(1:3, b$pe, pch = 16)
lines(1:3, b$pe)
segments(x0 = 1:3, x1 = 1:3, y0 = b$lb, y1 = b$ub)

#axis2
axis(2, las = 1, cex.axis = 0.8 , tck = -0.005)
mtext(side=2, text="Pr(Negotiation)", line = 1.25 , cex = 0.7)

#axis1
par(mgp = c(0,-0.5,0))
axis(1, at = c(1,2,3), labels = c("Weaker/ \n Much Weaker", "At Parity", "Stronger/ \n Much Stronger"), cex.axis = 0.7 , tck = -0.005, padj = 1)


text("Bilateral", x = 2.5, y = 0.075, cex = 0.8)
text("Mediated", x = 2.5, y = 0.3, col = "salmon", cex = 0.8)
dev.off()

#Figure 7: Distances ####
rm(list = ls())

## main_data ####
d <- readxl::read_excel("data_application.xlsx")

## model for simple application  ####
d$dist2 <- d$dist^2

d$int1 <- d$dist*d$territorial
d$int2 <- d$dist2*d$territorial

m1   <-  glm(formula = PeaceAgreement ~ dist + dist2 + 
               Intensity_rate + 
               RebStrBin.imp + 
               duration_dy_ln + imp.pop.ln + imp.rgdppc.ln +
               v2x_polyarchy + rebno_ct +
               territorial + int1 + int2, 
             data = d,
             family = binomial(link = "logit"))

summary(m1)

# coefficient draws for simulated CIs
coef.draws <- mvrnorm(n = 10000, mu = m1$coefficients, Sigma = vcovHC(m1, "HC1"))


### distance 
kms <- seq(0, 10000, by = 100)

## model frame
mf <- model.frame(m1)

## model matrix
mx <- model.matrix(m1)


# feed data frames
#colnames(coef.draws)
a <- 0
feedin1 <- cbind(Intercept = 1,
                 dist  = kms, 
                 dist2 = kms^2,
                 Intensity_rate = median(mf$Intensity_rate),
                 RebStrBin.imp = median(mf$RebStrBin.imp),
                 duration_dy_ln = median(mf$duration_dy_ln),
                 imp.pop.ln = median(mf$imp.pop.ln),
                 imp.rgdppc.ln= median(mf$imp.rgdppc.ln),
                 v2x_polyarchy= median(mf$v2x_polyarchy),
                 rebno_ct = median(mf$rebno_ct),
                 territorial = a,
                 int1 = a * kms,
                 int2 = a * kms^2)
a <- 1
feedin2 <- cbind(Intercept = 1,
                 dist  = kms, 
                 dist2 = kms^2,
                 Intensity_rate = median(mf$Intensity_rate),
                 RebStrBin.imp = median(mf$RebStrBin.imp),
                 duration_dy_ln = median(mf$duration_dy_ln),
                 imp.pop.ln = median(mf$imp.pop.ln),
                 imp.rgdppc.ln= median(mf$imp.rgdppc.ln),
                 v2x_polyarchy= median(mf$v2x_polyarchy),
                 rebno_ct = median(mf$rebno_ct),
                 territorial = a,
                 int1 = a * kms,
                 int2 = a * kms^2)


# dimension checks
dim(coef.draws)
dim(t(as.matrix(feedin1)))

# ystar for first plot
ystar.1 <- coef.draws %*% t(as.matrix(feedin1))
range(ystar.1)
yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
range(yhat.sim1)

# point estimates and cis
yhat.sim1.pe <- apply(yhat.sim1, 2, quantile, probs =  0.50)
yhat.sim1.lb <- apply(yhat.sim1, 2, quantile, probs =  0.05)
yhat.sim1.ub <- apply(yhat.sim1, 2, quantile, probs =  0.95)

## color ####
percent <- 60
fark <- 60
col.ref <- rgb(red = col2rgb("gray80")[1] - fark, 
               green = col2rgb("gray80")[1] - fark/2,
               blue = col2rgb("gray80")[1], max = 255,
               alpha = (100-percent)*255/100)

## plot 1 ####
dev.off()
png("distance.png", units="cm", width = 16, height = 14, res = 600)
par(mfrow = c(1,2), mar = c(2.5,2.5,3,1), oma = c(0, 0, 0, 0), mgp = c(1.5,0.15,0))
plot(kms, 
     yhat.sim1.pe, t = "l", 
     ylim = c(0, 0.35), col = "black", bty = "n",
     axes = F,
     main = "",
     xlab = "Distance to Capital (km)",
     ylab = "Pr(Agreement | Negotiation = 1)",
     cex.lab = 0.75,
     cex.main = 0.75)
lines(kms,  yhat.sim1.pe, t = "l", lwd = 2)
lines(kms,  yhat.sim1.lb, t = "l", lty = 2)
lines(kms,  yhat.sim1.ub, t = "l", lty = 2)
polygon(x = c(kms, rev(kms)),
        y = c(yhat.sim1.lb, rev(yhat.sim1.ub)), 
        col = col.ref, border = NA)
axis(2 , las = 1, cex.axis = 0.6, tck = -0.01)
par(mgp = c(0,0.1,0))
axis(1, cex.axis = 0.75, tck = -0.01)
mtext("Governmental Conflict", outer = F, cex = 0.75)

## plot 2 ####

# ystars
ystar.2 <- coef.draws %*% t(as.matrix(feedin2))
range(ystar.2)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)
range(yhat.sim2)

yhat.sim2.pe <- apply(yhat.sim2, 2, quantile, probs =  0.50)
yhat.sim2.lb <- apply(yhat.sim2, 2, quantile, probs =  0.05)
yhat.sim2.ub <- apply(yhat.sim2, 2, quantile, probs =  0.95)

par(mgp = c(1.5,0.15,0))
plot(kms, 
     yhat.sim2.pe, t = "l", 
     ylim = c(0, 0.35), col = "black", bty = "n",
     axes = F,
     main = "",
     xlab = "Distance to Capital (km)",
     ylab = "Pr(Agreement | Negotiation = 1)",
     cex.lab = 0.75,
     cex.main  = 0.75 )
lines(kms,  yhat.sim2.pe, t = "l", lwd = 2)
lines(kms,  yhat.sim2.lb, t = "l", lty = 2)
lines(kms,  yhat.sim2.ub, t = "l", lty = 2)
polygon(x = c(kms, rev(kms)),
        y = c(yhat.sim2.lb, rev(yhat.sim2.ub)), 
        col = col.ref, border = NA)
axis(2 , las = 1 , cex.axis = 0.6, tck = -0.01)
par(mgp = c(0,0.1,0))
axis(1, cex.axis = 0.75, tck = -0.01)
mtext("Territorial Conflict", outer = F , cex = 0.75)
mtext("Negotiation Venue and Peace Agreements", outer = T , cex = 0.85, line = -1.5, font = 2)
dev.off()

